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Abstract 

We propose a generalization of the cavity method to quantum spin glasses on fixed connectivity 
lattices. Our work is motivated by the recent refinements of the classical technique and its potential 
application to quantum computational problems. We numerically solve for the phase structure of a 
connectivity q = 3 transverse field Ising model on a Bethe lattice with it J couplings, and investigate 
the distribution of various classical and quantum observables. 
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I. INTRODUCTION 



The appearance of finite connectivity trees {eg Cayley trees, Bethe lattices) in the study 
of spin glasses has a long history that goes back to the first papers on the Sherrington- 
Kirkpatrick modelii^. Thouless, Anderson and Palmer- showed by means of a diagram- 
matic expansion of the partition function that in the spin-glass phase the mean-field the- 
ory is defined by the diagrams which describe an infinite tree with connectivity q ^ 1. 
They then wrote the mean field equations^^for such a tree and simplified the results in the 
small coupling, large connectivity regime (since in the SK model q = N and the couplings 
Jij oc 1/V^)- Those mean field equations are known as the TAP equations and the pecu- 
liar characteristics of their solutions, in particular their large number-, were an important 
indicator of the complexity of the spin glass phase that arose in parallel with the remarkable 
developments in the study of the replicated free energy^i^. 

Many authors have since studied the problem of a spin glass on the Bethe lattice^'^'^'i^i'ii 
with two distinct motivations. The first of these has been to attempt to find a model of a 
short ranged spin glass where one can rigorously assess Parisi's picture for the organization 
of the Gibbs states in short-range systems^. The second has come from computer science 
whereby a set of optimization problems can be recast as frustrated problems on random 
graphs with the local connectivity of a tree^^. 

Despite much work, the Bethe lattice has not yielded a decisive verdict for or against the 
Parisi picture of a multitude of Gibbs states in the ordered phase of a spin glass although the 
case for it is perhaps stronger here than on regular lattices with short ranged interactions. 
This is tied up with the question of defining the infinite Bethe lattice limit starting from finite 
graphs. It is possible to do so either via a sequence of Cayley trees with random, frustrating, 
boundary conditions or via a sequence of random graphs with fixed connectivity. The latter 
sequence has frustrating loops of typical size diverging 0(log(A^) / log(g — 1)) for a graph with 

points and connectivity q and thus locally looks like a tree^. The Cayley tree sequence 
can always be analyzed in terms of a recursion relation that we review below and does not 
appear to lead to the Parisi structure^. However, it has been argued that the random graph 
problem is fundamentally different and does lead to replica symmetry breakingi^. 

Happily, the perspective provided by spin glass theory on optimization problems has been 
quite fruitful regardless. Starting in the early SO's^*^, it became clear that much was gained 
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by the recognition that various optimization problems in computer science were equivalent 
to finding ground states of certain statistical mechanics problems. A typical example of 
this connection is given by the k-SAT problem, which asks "Given a boolean expression J 
on bits cTj composed of the conjunction of M clauses, each of which involves exactly k 
of the bits, is it satisfiable by some bit assignment?" This can be recast in Hamiltonian 
form by writing a cost function iyj[{(Tj}] that evaluates the number of violated clauses. In 
this language, the bits naturally become Ising spins, the expression J becomes a particular 
instance of some spin glass and the original question requires determining the ground state 
energy. The large A^ limit is a problem in statistical mechanics and typical-case analysis 
for k-SAT becomes the disorder averaged analysis of spin glass theory. In this limit, /c-SAT 
develops several phase transitions as a function of a = M/N, the number of clauses per 
hiir^. With increasing a, the most salient features are that the problem goes from being 
easily solved and satisfiable, to an intermediate glassy phase with many local ground states, 
to a typically unsatisfiable phase where the ground state energy density is positive. 

A key role in these developments has been played by the so called cavity method — a 
complex of analytical and numerical techniques refined recently^^^^^ — for studying classical 
spin glasses on tree-like graphs. Applied to the k-SAT problem, the cavity method suggests 
the above phase diagram and provides numerical values for its critical points. Moreover, the 
technique can be applied to a particular instance of k-SAT and the information so obtained 
about the free energy landscape now guides the search procedure in state-of-the-art k-SAT 
algorithms'^. 

In this paper we turn to quantum spin glasses on Bethe lattices, specifically to the problem 
of extending the cavity method to their analysis. As in the classical case, there are two 
distinct reasons to be interested in these systems. There is the intrinsic interest of the 
interplay between quantum mechanics and spin glass behavior, about which the difficulties of 
the classical case serve as both caution and enticement. A second motivation now arises from 
the rapid recent developments in quantum computing. In particular, the discovery of the 
ground state of a classical spin system Hj derived from a computational problem is ideally 
suited to solution by the adiabatic algorithm^, which allows a quantum computer to solve 
such a problem by means of an adiabatic change of the parameters in the Hamiltonian. The 
adiabatic algorithm would trace a path in operator space starting from a simpler Hamiltonian 

and ending at Hj. For instance, Hq could be the Hamiltonian for independent spins 
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in a large transverse magnetic field Bt and the path could be to slowly lower the field and 
raise the Ising couplings J. Starting from the easily found ground state of Hq a sufficiently 
slow evolution will carry the system to the ground state of Hj. The power of this algorithm 
is therefore measured in the scaling of the evolution time with the number of variables A^. 

This discussion certainly suggests that a careful study of the phase transitions encountered 
on the path between Hq and Hj is in order. However, the necessity of understanding the 
"deep" quantum spin glass phase far from the phase transition is also clear. The structure 
of this phase, especially its nontrivial energy landscape, is probably as important as the 
nature of the phase transition. In this paper we take a first step toward understanding 
the quantum spin glass phase by generalizing the cavity method used to study the classical 
problem. Compared to the replica method or direct study of the quantum TAP equations, we 
believe the cavity method outlined in this paper provides much more physically transparent 
information about the SG phase. It also may be applied to many other quantum phase 
transitions on the Bethe lattice, whether induced by disorder or more typical symmetry 
breaking. 

We note that the topic of quantum spin glasses on trees has been tackled before in 
the literature from the point of view of statistical mechanics^ and computation theory^!. 
However, in order to proceed analytically and due to the complexity of the problem, these 
works have employed several uncontrolled approximations. We will compare with them 
briefly in the core of the paper. While we have used approximations in our (numerical) 
analysis as well, we believe they are better controlled, since we show how to systematically 
improve them and how the results are robust with respect to the improvements. More 
evolved computational methods for tackling very similar problems have appeared in the 
literature on dynamical mean field theory— and could possibly be applied to the quantum 
cavity method on glasses. These methods may allow a study of the zero temperature case, 
which our approximations could not capture. 

Let us turn now to a very brief overview of the results in this paper. We first formulate 
the cavity method for the transverse Ising spin glass on a Bethe lattice. This is a Markov 
process for the on-site effective action whose stationary probability distribution can be found 
numerically using a population dynamics algorithm. This will be the main result of the 
paper. In this way we obtain the phase boundary in the {Bt,T) plane, the values of the 
usual thermodynamic quantities (free energy, energy, entropy and Qea, the Edward- Anderson 
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order parameter—) and more typical 'quantum' quantities like the single spin von Neumann 
entropy. Indeed, from the fixed point probability distributions for the effective actions we 
can calculate distributions for all of the statistical properties of the system. This leaves 
much room for further work. 

The plan of the paper is the following: we introduce the classical cavity method for 
uniform ferromagnetic systems and spin glasses in Section [Tll we propose its generalization 
to the quantum spin glass problem in llllj we apply it to the numerical study of a Bethe 
lattice with connectivity three in IIV[ Discussion and directions for further work will be 
presented in the last sections. 

II. CLASSICAL CAVITY METHOD 

The cavity method is a way of finding the spin glass free energy by means of self- 
consistency equations for the probability distributions of quantities characterizing the statis- 
tics of the spins (the cavity fields) . It has several virtues compared with the replica method, 
in particular if applied to spin glasses with finite connectivity. We will return to them after 
we have explained how the cavity method works. 

A. The Bethe-Peierls method for the ferromagnet 

For conceptual clarity we will first recall the "cavity method" in its simplest form as 
applied to a ferromagnetic Ising model living on a Bethe lattice with connectivity q. As 
mentioned above, there is considerable subtlety in defining a sensible infinite tree model, 
but for this simple case, a sequence of Cayley trees with uniform small boundary fields to 
break the Ising symmetry will suffice. In any event, we shall focus here on the local structure 
of the model. 

The classical Hamiltonian is 



where cxj G {±1} are Ising spins, Jij = J > is a uniform ferromagnetic coupling and the 
sum is over bonds in the Bethe lattice. To create a cavity, we pick a spin ctq in the Bethe 




(1) 
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lattice and imagine removing it: 





-1^ 



(2) 



Each of (To's neighbors is a cavity spin, connected to g — 1 spins and sitting at the root of 
a branch of the original tree. Notice that in the absence of ctq these q branches are entirely 
independent. Similarly, a cavity spin ai mediates the only interaction between the g — 1 
branches sitting above it. 

We can define three important operations we can perform on graphs with cavity spins: 
iteration, merging, and link addition (Fig. [H) For our immediate purposes, the most 
important of these is iteration, which takes q — l rooted branches of depth / and links them 
into a single new branch of depth / + 1 with a new cavity spin at the root. Thus we can 
construct an arbitrarily large Cayley tree by iterating inward from its boundary spins many 
times followed by a final merge operation to form the center. 



FIG. 1: The three cavity operations. Although we have not labeled them in the Figure, the new 
links are Jqi, Jo2i etc ... 

Consider the iteration operation: the added spin ctq receives thermodynamic information 
regarding each of the q — 1 branches only through the thermal distribution of the cavity spins 
(Ti, (Tq-i. In the absence of ctq, each of these Ising variables has independent statistics 
characterized fully by its thermal probability distribution: 



which defines the cavity field hi. Since ai G {±1}, there are only two possible configurations 
of the spin and only 2 — 1 = 1 real numbers are needed to characterize the probability 






(3) 
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distribution. In this sense, the cavity field hi is merely a good parameterization of the 
distribution ipi^ai). We emphasize this viewpoint because it will naturally generalize to the 
quantum case. 

We introduce a simple graphical convention for cavity spins: 

(4) 

The open circle indicates a spin variable and the wiggly line indicates the effective field 
attached to it. With this notation, an iteration operation can be represented: 

(5) 

where the filled circle indicates summing out a spin variable. More formally, the state of the 
spin (Jo depends on the state of the g — 1 spins as 



^o(cro) = ^ X] exp I /3 ^ JoiO-Qcr* j 

o-l,...,o-„_i=±l \ i / 



i/^i(ai)...V^5_i(ag_i), (6) 



where Z is a normalization factor so ^^^^ ^(o"o) = 1. In terms of cavity fields this equation 
implies 

1 

hQ = ^Y^ tanh-\tanh(/3Joi) tanh(/3/i,)) = U{{hi], {Jo,}). (7) 
^ i=i 

Finally, to solve the ferromagnet, we note that all the Jij = J and apply uniform boundary 
fields (say slightly positive) to the Cayley tree "leaves". From this uniform starting point, 
we expect to find fixed points for the cavity fields under iteration given by 

h = tanh-i(tanh(/3J) tanh(/?/i)) (8) 

P 

This is precisely the Bethe-Peierls self-consistency equation for a mean field ferromagnetic 
in a lattice of coordination number q. 

B. Classical Spin Glass on a Cayley Tree 

The careful reader will have noticed that until Eq. ([8]), we did not anywhere exploit the 
uniformity of Jij or hi in the foregoing discussion. With this foundation laid, we can make 



short work of the classical Bethe lattice Ising spin glass. Again, we consider a Bethe lattice 
defined as a limit of Cayley trees with fixed boundary conditions rather than as a limit of 
random graphs. The Hamiltonian is now 

Hj = -Y^ Jijaiaj (9) 

where the Jij are i.i.d. random variables drawn from some distribution P{J)- For simplicity 
and because of its connection to computational problems, we will restrict our attention to 
the ±J model 

although much of the discussion has broader validity. 

The iteration equation ([6]) is still valid for particular realizations of the J, but since these 
are random variables, it now defines a Markov process for the cavity fields. Throughout the 
graph, these fields will be site dependent random variables, but deep inside the tree, they 
ought to be distributed according to a probability distribution P{h) that represents a fixed 
point of the Markov process.— This fixed point distribution will satisfy 

W dKP{h,) {6{h - U{{h,}, {Joi}))j . (11) 
1=1 

In terms of the spin distribution, this becomes the functional equation 

l[Di^,p[4^,] (sma) - M<^; m, M)]),; (12) 

where ipoi^') is given by Eq. (j6]). 

This distribution is the order parameter for the spin glass. It is a 5 function at /i = in 
the high temperature phase. As the temperature is lowered, P{h) broadens to have finite 
support below some mean field-like phase transition. Defining r = T/J, this phase transition 
is located at = 1/ tanh~^(-7==). It is possible to write the free energy per site in terms 
of P{h) as 

/q p 2 

n dh,P{h,)Fg+i - ^ / n dh^P{hi)F2, (13) 
i=i ^ i=\ 

1 / ^ ^ \ 

Fg+i = -- Mn ^ exp(/?^ JoiO-QO-i) JJ?/'i(o-,) \ , (14) 

= (in ^ exp(/5Ji2ai 0-2)^^1 (cri)V^2(o-2) y • (15) 
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This rather comphcated looking expression is actually simply the average change in free 
energy due to a merge operation minus | times the average change in free energy due to a 
link addition. Graphically, 



F ={ T 



2 



(16) 



where JF of a diagram is the free energy of a system with spin variables given by the unfilled 
circles. 

It is possible to see that if P satisfies (ITT!) then 6F/6P{h) = 0. Other expressions for 
the free energy F have appeared in the literature but it has been shown that they are all 
equivalent to each other, if the consistency equation (ITTl) is verified. The role of the free 
energy in the cavity method is secondary as one does not solve the variational problem, as in 
the replica method, by working on the free energy directly. Rather, one finds the probability 
distribution P by analytical or numerical methods and then derives all of the statistical 
observables from P.— The equivalence of the two formulations has been put forward in Ref. 
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151 and in many other works. 
We have up to this point assumed that the underlying lattice is in fact a tree and that 
the removal of a spin to create a cavity completely decouples the neighboring branches. On 
such models, the cavity method we have described is exact. An important generalization of 
the cavity method arises in its heuristic application to locally tree-like random graphs where 
the typical loop size diverges logarithmically and any finite neighborhood is a simple tree. In 
this case however, the decoupling of cavity spins is not exact and it is necessary to introduce 
the so-called "replica symmetry breaking" ansatz on the structure of thermodynamic states 
in order to correctly describe the frozen (glassy) phase. Although this is a vital component of 
the modern understanding of such glasses, we have not included it in our quantum treatment 
below because we believe their is still much to be understood about the quantum model 



15 for a more detailed 



on the much simpler trees. We refer the interested reader to Ref. 
description of replica symmetry breaking in the classical cavity method. 

Now, having introduced the cavity method in classical statistical mechanics we go on to 
generalize it to quantum mechanics. 
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III. QUANTUM CAVITY METHOD 
A. Exact Framework 

We consider the modification of tlie Hamiltonian ([T]) due to tlie introduction of a trans- 
verse magnetic field 



This is called the transverse field Ising spin glass in the literature. The Ising variables of the 
previous section have been replaced by Pauli matrices and the magnetic field couples to 
the matrices a^. The fact that and do not commute gives rise to a host of interesting 
new features due to the interplay of quantum mechanics and disorder—. 

The usual Suzuki- Trotter decomposition allows us to rewrite the problem in terms of Nt 
Ising spins per quantum spin, where the number Nf needs to be sent to infinity eventually. 
The additional dimension which is introduced in this way is the usual imaginary time. 
The CTiCTj interactions are time-translation invariant (the disorder is correlated in the time 
direction) while the terms give a ferromagnetic nearest-neighbor interaction in the time 
direction. Before writing the Hamiltonian let us introduce some notation. 

For any finite Nt we will refer to the Ising spin configuration at a given site z as a "rod" 
of spins. The rod at site i is described by A^^^ spins cri(t) where t takes values from to /3 
in steps of At = (3/Nt, with periodic boundary conditions cr(0) = (r{(3). This notation is 
convenient if the limit Nt ^ oo is eventually performed, since the rod is represented by a 
function cr(t) : [0,P] — > { — 1, 1} with crj(O) = <Ji{P). The rod statistics are described by a 
probability distribution ?/;[(T(t)], a functional of o"(t), that gives a positive real number for 
every configuration a(t). The normalization condition reads X]{o-(t)} '^['^(^)] — ^■ 

The partition function is written as 




(17) 




(18) 



where the Hamiltonian is: 




(19) 



t (ij) t i 

We can also write this as a sum over links of the energy of a link: 




(20) 
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where a fraction 1 /g of the imaginary-time interaction is associated to each hnk (there are 
q hnks per spin). Here ^ = \ \ncoi\i{(3Bt/ Nt). 

Notice that F > so the system is ferromagnetic, and moreover when jSBt/Nt ^ 1 then 
r ^ 1 and the couphng along the time direction is strongly ferromagnetic. In particular, 
for Bt = the spins in any given rod are locked together as Nt useless copies of a single 
Ising spin. Thus, the results reduce to the classical case smoothly. 

The spatial tree-like structure of the original problem is reflected in the tree-like structure 
of the interaction between rods. We can therefore imagine an iteration process with rods 
replacing the spins, in which we have q — 1 cavity rods ipilaiit)] which are merged and 
determine the state of the rod ^/'o[c"o('f)] (see Fig. [2]). This corresponds to a recursion 
relation for the calculation of the partition function of the branches, analogous to the classical 
equation IQ: 

MMt)] = ^ E e-s-^^«n^4^.W]- (21) 




FIG. 2: Iteration of cavity rods. There are periodic boundary conditions in the imaginary time 
(vertical) direction. Other cavity operations are analogous. 

Just as in the classical case, this iteration equation is already enough to solve an inter- 
esting ferromagnetic problem. Consider the case q = 2 and Jij = J. For q = 2 our Bethe 
lattice is a simple chain of spins, moreover all J's being equal we can take J > without loss 
of generality. We then recover the well known ferromagnetic Ising chain with a transverse 
field, an exactly solvable system. One way to solve it is to use Onsager's transfer matrix 
method. In fact, the iteration equation ( l2Ti) can be rewritten as just such a transfer matrix 
equation where ip is the 2^ dimensional vector and e~^^°'^ the transfer matrix T: 

^^o = T-^i. (22) 
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The fixed point of tliis iteration gives the eigenvector ip corresponding to the largest eigen- 
value Z. In the limit Nt — > oo, this is an exact solution which contains all the information 
about classical and quantum phase transitions. 

Returning to the case of the spin glass, we can write down the quantum cavity fixed point 
equation analogous to the classical equation fllip immediately: 

(23) 

= / (rc=i D^.Pfp [^P,] rfJo.P( Jo.)) 6 [^[a(t)] - 7/>o[a(t); {Jo., 

where the iterated rod action ipQ[a{t); {Joi, ipi}iZl] is given by Eq. fl2T]) and P(Joi) is the fixed 
prior distribution for couplings Eq. (ITUIl . This is a functional equation for Ppp [ip[a(t)]], the 
fixed point probability distribution of the effective distribution describing iterated cavity 
rods. In the limit A^t ^ cxd, it is exact but difficult to solve in closed form. It is certainly 
possible that analytic progress can be made, but we have not succeeded thus far. However, it 
is amenable to numerical study at finite Nf under certain approximations and also perhaps 
by continuous time Monte Carlo for Nt ^ oo. In the remainder of this section we will 
explore the finite Nt approach. 

We must first parameterize our generic vector ip in the 2^*-dimensional space of the 
configurations of the rods. In principle it is described by 2^* — 1 real numbers, which can 
be reduced by a factor 0{Nt) by exploiting time-translation symmetry and the periodic 
boundary conditions. A natural way to parameterize it is in term of the effective action of 
the rod 

^[a(t)] = e-^["l (24) 
where we expand S' in a series of increasing clusters of interacting spins. 

S[a]= - log Z -hAtJ2(^it) -J2 ^^'C'^') {t' - t)a{t)a{t') - 

t t,t' 

- J] At3c(3)(t' _ t,t" - t')a{t)a{t')a{t") + ... (25) 

In principle, the sum includes up to A^^-spin interaction terms (the normalization factor has 
been included as a spin-independent term in the effective action). In practice, we truncate 
the action expansion at second order to keep the numerical requirements manageable. We 
comment below on the limits in which this truncation is exact. Notice that C*^^-'(t) is the 
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kernel for 2-point in time interactions, not the dynamical two-point correlation function, 
often denoted c^^\t) = {a{t)a{0)). 

The functions h, C'-*-' are random quantities characterized by the Markov process defined 
by the iteration procedure. By writing the representations of the vectors tp in terms of the 
effective action fl25p we can rewrite the iteration equation as 

q-l 

^-S[a,{ho,Co}] ^ ^ g-E^Ji /3^^o, J~|'g-5[a,{/.„C,}]_ (^26) 

{al{t)},...,{a,-l(t)} j=l 

(2) (S) 

This gives an implicit update map from the 'old' q — 1 parameters hj, C- , Cj , ... and the 

(2) (3) 

couplings Joj,Bt to the 'new' parameters /107C0 7 Cq > •••• The statistics generated by this 
Markov process and in particular its fixed point distribution 

Pih,C^^\C^'\...) (27) 

are the solution of the problem. 



B. Approximations 

1. Action Representation 

By parameterizing the quantum dynamics through the action expansion (125!) we do not 
make any a priori assumptions about the nature of the spin-spin correlations in the time 
direction - that is, if we could keep all of the terms in the expansion, it would be an exact 
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treatment. We avoid in this way the spherical approximation which has been used in Ref. 
since the results thus obtained do not reduce to the well-known classical results for Bf = 0. 

In practice, however, we truncate the cavity actions to second order. Most usefully, this 
corresponds to the leading order term in a large connectivity expansion of the effective 
action. Indeed, in a large q treatment, in which the couplings Jij must be scaled as \l ^fq 
for the disordered model, one finds that the one- and two-body terms in the effective action 
for the rod at the root of a tree are O the three- and four-body terms are O (l/g^) 

and so on^. The truncation to second order is thus both exact and necessary in the q ^ 00 
limit. We note that the oft-used static approximation^ is thus incorrect even at infinite 
connectivity for disordered models. The truncation to second order is also exact at high 
temperature, regardless of the value of q. 
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Numerical investigations of small systems suggest that the higher order interactions are 
quantitatively small more generally, even at g = 3. This is especially true at small Bt, 
where the strong ferromagnetism in the time-like nearest neighbor bonds dominates. We 
note that there is also some error in the numerical fitting of an iterated cavity distribution 
to a truncated action. This problem is reminiscent of a maximum entropy model for the 
statistics of signals and the inverse Ising problem of computer science^I. However, even in 
the disordered system, the interaction of the spins along imaginary time is always ferromag- 
netic and the problem does not present the difficulties that usually accompany the inverse 
problems in general statistical mechanics. In other words, although the original problem is 
fundamentally frustrated, none of this frustration appears in the single spin dynamics. The 
frustration is taken care of in the treatment of the relevant parameters C*-*^ as random 
numbers. 

For comparison, we note that our approach is closely related to self-consistent dynamical 
mean field theory (DMFT) methods, which also truncate to non-local two point interactions 
in the single-spin effective action (these are the so-called Weiss functions or bare Green's 
functions of DMFT treatments). However, DMFT techniques cannot be used straightfor- 
wardly for disordered systems. One cannot assume that the imaginary-time propagator (or 
action) for a spin is equal to that of the spins which surround it but rather that they come 
from the same probability distribution. If this observation is taken into account then we 
expect to recover functional equations for the distribution of single site Green's functions 
analogous to those we have proposed for the single site action. 

2. Finite Discretization Error 

Finally, we comment on where in the (r, Bt) phase diagram the fixed Nt approach will 
give reliable results. As mentioned above, on the classical line i?t = 0, all spins in a rod 
are locked together by divergent nearest neighbor interactions and the Trotter decomposed 
system (for any Nt) reduces to the classical system exactly. In the opposite limit of large Bt 
at fixed Nti the planes of the Trotter system decouple as F ^ 0. Each of these planes is an 
exact copy of the original Bt = spin glass model except with couplings Jij/Nt. These will 
therefore undergo independent thermal phase transitions at r = Tc/Nt and no r = critical 
field will be detected. This explains the phenomenon of "asymptotic critical lines" that we 
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note in our finite Nt phase diagrams. 

More generally, for the finite Nt approximation to be valid, Nt > ktyp where fctyp is the 
typical number of kinks in a rod. For a single spin in a transverse field, a straightforward 
calculation shows (k) = 2j3Bttanh. /SBt, which for large /3Bt reduces to (k) ^ 2(3Bt. Within 
the paramagnetic phase this calculation remains nearly exact, although in the spin glass 
phase {k) will be suppressed by the presence of longitudinal cavity fields. However, this 
inequality A^^ > 2f3Bt remains a good indicator of the quality of the approximation and 
agrees with the regions where it clearly breaks down. 

Thus, this expansion is particularly useful close to line -B^ = 0, and we will see that 
it gives reasonable results and insights into the structure of the problem also deep in the 
quantum spin glass phase. On the other hand, the other region of interest T = Bt — B^^^^ 
should not be addressed with this expansion, unless the description of the spins in terms of 
continuous time functions a{t) turns out to be treatable in the future. In this paper, we will 
not be able to make definitive statements about the nature of the quantum phase transition 
which occurs at this point but we are definitely able to make statements about the nature 
of the spin glass phase when quantum effects are not negligible. 

In the next section we show how a simple minded Trotter discretization with Nt relatively 
small {Nt = 6 to 11) delivers a great deal of information about the quantum spin glass phase. 

IV. NUMERICAL RESULTS 

We solve the fixed point equation (123|1 numerically using a population dynamics algorithm 
analogous to that of Mezard and Parisi^. We represent P[ip] by a finite population of Nrods 
rod actions {^/'i[(T(t)] = e"'^''"'^*-"}, where the expansion fl25|) of Si is truncated to second 
order. Each rod action is therefore specified by 1 + [^J distinct numbers {h,C^'^\t)) after 
exploiting the periodicity of imaginary time. 

This population is initialized from an appropriate uniform distribution and then iterated 
as follows: 

1. Select q — 1 rods ipi randomly from the population and q — 1 random Joj. 

2. Use ( l2Ti) to calculate the effective action on an iterated cavity spin ipo from these 
rods. In principle, higher order interactions may be generated but we truncate them 
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by finding the second order action that exactly reproduces the free energies of a series 
of domain wall configurations of varying width. 

3. Randomly replace one element of the population with ipQ. 

4. Repeat until convergence in some measure of the population, for example the order 
parameter qea- 

In practice, this procedure converges quickly deep in either the glassy or paramagnetic phase 
but slows near the phase transition as the flow of P[ip] under the iteration equation near 
the paramagnetic fixed point becomes marginal. 

Given Ppp[ip], we can calculate the sample averaged free energy density and local ob- 
servables such as the link energy, site magnetization, Edwards-Anderson order parameter 
Qea = ((^Ot)^ reduced von Neumann entropy S^onN = — trpo log2(po) by standard 
Monte Carlo sampling of these quantities. The free energy density is given by equation f|T6l) . 
The reduced entropy and transverse magnetization may be derived in the usual way from 
the reduced density matrix po for a spin. We calculate this by performing a merge operation 
onto a "broken rod" (see Fig. [3]), in which periodic boundary conditions are not enforced. 
The various elements of the reduced density matrix correspond to imposing different values 
on the top and bottom spin of the broken rod and summing out all other spins in the parti- 
tion function. Finally, calculating the average internal energy can be done by averaging the 
Hamiltonian ( |T7I) . 

A. Numerical Results 

1. Phase Diagram 

We present numerical results for an investigation of the q = 3 connectivity model using 
a naive (exact) approach to the exponential summation involved in the cavity iteration and 
merging operations.— Fig. Hl^a) shows the phase diagram calculated at Nt = 10, Nrods = 
2500, Niter ~ iOOONrods and suggestively fit to A''^ — oo using asymptotic expansions in 
l/N^. Qualitatively, all is as might be expected: 

• At any Nt, the phase transition curve predicts a Bt = critical temperature in agree- 
ment with the analytic prediction of = 1/ tauh~^{l/y/q — 1) ^ 1.13. 
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FIG. 3: Merging of two cavity rods onto a "broken" central rod. Here A'^t = 3 but instead of 
imposing periodic boundary conditions everywhere we depict cr(0) and cr{Nt) as independent spins. 
Each vertical link corresponds to a F coupling while each horizontal link corresponds to a Jij/2Nt 
coupling. The dashed lines connecting the cavity rods indicate identification of the top and bottom 
spins for each rod and the wiggly line indicates an effective action. By imposing different values on 
the top and bottom spin of the broken rod and summing out the rest, the elements of the reduced 
density matrix pQ may be determined. 



The upturn in the Nt = 10 phase boundary at low temperature is due to the finite 



discretization of time, which leads to an asymptotic phase transition line at r = - 
0.113. 

• While the fits to A'^t ^ oo are certainly approximate, we believe that the true r = 
critical field lies between 1.5 and 2. We believe continuous time techniques will allow 
dramatic refinement of this estimate and investigation of the quantum critical region. 

Our phase diagram clearly disagrees with that of Kopec and Usadel^S^ who treat the identical 
model using a soft spherical approximation and find that both the critical temperature and 
critical transverse field are depressed relative to our values. Presumably this suppression of 
ordering arises due to the stronger effect of fluctuations in the softened model. 

Figure 111(b) shows the instance averaged single site von Neumann entropy S^onN which has 
a remarkably clear maximum near the phase transition curve above the classical line. This 



reflects t 
See Ref. 



le strength of quantum correlations even at the finite temperature phase transition. 
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for discussion of local measures of entanglement at finite temperature. 
Zooming in on the horizontal stripe at Bt = 1 indicated on the phase diagram, we 
find that qea vanishes linearly at the critical temperature (Fig. [5t^a)). This reflects the 
underlying broadening transition in P[ip], which can be seen sharply in the variances of 
each of the effective action coefficients (Figs. [5](b,c)). We use this behaviour to estimate 



17 



sharp transition points despite softening due to critical slowing in the convergence of our 
procedure. 

Finally, we note that much of the phase diagram is surprisingly stable to variation in Nt. 
We have explored various regions of the phase space at A^^^ = 6, 7, 8, 9, 10, 11. The classical 
line {Bf = 0) at all temperatures is completely stable down to Nf = 1 as expected. Perhaps 
more surprisingly, moving between Nt = 8 and Nt = 10, qea is essentially stable below 
Bt = 1 down to temperatures r ~ 0.15. Of course, the high field, low temperature part 
of the phase transition curve moves downward as the finite discretization asymptote goes 
towards the r = axis. See Fig. [6] for the low temperature critical curves estimated using 
vertical stripes run at five different temperatures (corresponding to P = 3.5,4,4.5,5,5.5) at 
various Nt. 
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FIG. 4: (a) Phase diagram at g = 3. The solid phase transition curve has been calculated at 
Nt = 10, Nj-octs = 2500, Niter = WOONj-ods on a fine mesh in the (r, Bt) plane. The vertical dotted 
line is the asymptotic critical line for large Bt at N = 10 {ie r = Tc/Nt). The points marked x 
with error bars indicate N ^ oo fits based on Fig. [6l The dashed transition curve is a weighted 
quadratic fit through the estimated low temperature points and the Nt = 10 points in the range 
0.5 < r < 1. This leads to an estimated Bf = 1.775 it 0.03. As this fit is clearly heuristic, we have 
suggested a much larger range for our estimate of Bf in the Figure. The stars and stripes indicate 
points in the phase space which we have investigated in more detail below, (b) The average von 
Neumann entropy SyonN (in bits) of a central spin as a function of (r, Bt) at N = 8. The dashed 
line indicates the estimated region of validity of the discretization approximation (Bt < NtT/2). 
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FIG. 5: These graphs ah correspond to the horizontal shce at i^t = 1 shown in Fig. [H The 
order parameter qea and all effective action variances undergo mean field like transitions at the 
critical temperature {eg. qea ~ I''" ~ ''"cl^)- This allows us to estimate the critical temperature 
precisely despite softening due to critical slowing down near the phase transition. These curves 
were calculated at Nt = 8. However, the results are stable to increasing Nt to 10 to within an error 
of ±0.005 in Tc. 
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FIG. 6: The phase transition curve in the high field, low temperature regime at various Nf. The 
vertical dashed lines indicate the finite discretization critical asymptotes. The estimated curve for 
Nt ^ CO in Fig. is given by fitting Bt{Nt) = a/N^ + b/N^ + c to the points calculated at each 
temperature and then sending Nt to infinity in the result. These fits suffer from a paucity of data, 
but are suggestive nonetheless. 
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2. Structure of the Glassy Phase 

Since the order parameter is a probability distribution of an action, there is a rich struc- 
ture to be investigated at even a single point (r, Bt) within the glassy phase. In Fig. [71 
the marginal probability distribution of the field term h in the cavity action is shown at 
the four points indicated on Fig. |4](a). The two lower distributions lie on the classical line 
{Bf = 0), one deep within the glassy phase and one near the transition. It is clear that the 
distinctive features of the classical solution are reproduced here: a Gaussian-like structure 
around h = near the phase transition with the appearance of delta function spikes on 
the integer fields deep within the phase. At Bt = I, the qualitative picture of spread from 
narrow Gaussian near the phase transition to broader, bumpier distribution remains. It is 
less clear whether the sharply defined spikes on integer fields would remain at r = with 
large Bt. 

Further structure can be found in the nontrivial probability distribution for the interac- 
tion terms that develop in the spin glass phase. Figure M shows the histogram for various 
marginal and /i-conditioned distributions of the nearest neighbor and next nearest neigh- 
bor interactions terms at (r = 0.25, Bt = 1) (cf. Fig. [Tiptop left)). We can qualitatively 
understand many features of these distributions: 

• The two-spin interactions are ferromagnetic and the effect of coupling to neighboring 
rods is only to enhance the ferromagnetic interaction from the bare nearest neighbor 
interaction on a single rod (F). Indeed, this F sets the minimum strength of C{6T), 
as can be seen in the top row. 

• The strength of two-spin interactions are strongly anticorrelated with the strength of 
the cavity field h as can be seen from the decomposition of the full marginal distri- 
butions of C{6t) and C{26t) into the small (middle column) and large (right column) 
cavity field conditioned distributions. Large cavity fields on a central spin come from 
large fields biasing neighboring rods. These fields pin the neighboring spins more 
strongly and reduce the ability of those spins to mediate interactions in time between 
the central rod spins, reducing the effective two-spin interaction. 

• The multimodal spikiness in these distributions reflect the spikiness in the low tem- 
perature cavity field distributions through the field-interaction correlation. 
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Unfortunately, we have yet to develop a more significant analytic understanding of these 
correlations nor a means to extrapolate them to zero temperature in the presence of the 
transverse field. 

Finally, we emphasize that the phase transition is signaled by a singular broadening of P 
rather than any singularity in its first moments or in the structure of the typical imaginary 
time action. This can be seen in the smooth evolution of (C(2At))j through the phase 
transition in Fig. [9^ and the similarly smooth evolution of the average two point correlation 
function {a{0)a{t)) j. This is in contrast to the ferromagnetic case in which the distribution 
P{ip) would exhibit spontaneous symmetry breaking but remain deterministic. 



P(h; T = 0.25, = 1 .0) P(h; x = 0.83, B, = 1 .0) 




-2-1012 -2-1012 



FIG. 7: The distribution of tlie field term of the cavity rod action at the four different points in 
phase space labeled by the stars on Fig. HI The distinctive distribution of the low temperature 
replica symmetric classical spin glass is reproduced in the bottom left corner while the three other 
points all lie closer to the phase boundary in the r or directions. 
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FIG. 8: Histograms of nearest neighbor (top row) and next nearest neighbor (bottom row) inter- 
actions in imaginary time at r = 0.25, Bt = 1 (in the SG phase). The second and third column 
provide the conditional distribution of the interactions given \h\ is small or large. 

V. DISCUSSION AND FURTHER WORK 

The cavity method has a long and illustrious history in the study of statistical systems - 
from Bethe's early work on the Ising ferromagnet to modern studies of random constraint 
satisfaction problems in computer science. In this paper, we have introduced a new variant 
for studying disordered quantum systems within an imaginary time formalism. This repre- 
sents an intuitively appealing, natural synthesis of the classical disordered model techniques 
with the quantum homogeneous models studied in DMFT. We note that our framework can 
be simply adapted to study many other transverse field Ising systems on trees with fixed or 
fluctuating connectivity - such as the ferromagnet, diluted ferromagnet or biased glass. 

We have shown that the transverse field Ising glass on a Bethe lattice of connectivity 
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FIG. 9: (a) The average next to nearest neighbor interaction in time C^'^\2At) of a cavity rod 
action at Bt = 1, varying the temperature. The bars indicate the variance of C(2)(2At). Notice the 
zero variance above the critical temperature, (b) Single site imaginary time correlation function 
{a{0)a{t)) J at T = 0.25 for various Bt passing through the transition at B^ ~ 1.85. 



three has a phase transition line all the way from the classical Bt = 0,T = 1.13 to the 
quantum Bt ~ 1.75, T = 0. At finite temperature, the transition is classical and mean field 
like. Inside the frozen phase the picture is similar to the classical case: when a randomly 
chosen spin is extracted from the graph, its effective action (analogous to the cavity field) is 
well-defined in the paramagnetic region but is a random functional in the spin-glass phase. 
In this phase all of the local observables, classical and quantum, are therefore also random 
variables. In principle one could also study the properties of the entanglement of distant 
spins, which cannot be done in the better known but fully connected SK model. 

The reader familiar with classical spin glass theory will have noticed our avoidance of 
the important issue of replica symmetry breaking (RSB), which is widely believed to be a 
feature of a correct treatment of the random graph Bethe lattice. In this connection, we note 
that our simpler treatment is indeed correct for models on Cayley trees with fixed bound- 
ary conditions. Quantum RSB phenomena certainly deserve further study: the conceptual 
difficulties of the RSB ansatz become even thornier in the presence of quantum tunneling. 
We note that formally breaking replica symmetry at the one step level (IRSB) should be 



straightforward in our framework. In analogy with Ref. 
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one should introduce popula- 



tions of populations of effective actions and weigh them according to their free energy (as 
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one does for different solutions of the TAP equations). This straightforward modification of 
the algorithm makes it computationally considerably more time consuming. For this first 
pass, we decided not to embark on such a project. 

Unfortunately, we are not aware of Quantum Monte Carlo or other numerical studies 
on the transverse Ising model on the Bethe lattice with which to compare our results. 
The primary difficulty that has prevented the direct simulation of this system is that it 
requires a very large number of spins to approximate the infinite system effectively. Since 
a random graph with fixed connectivity has an extensive number of loops of lengths i ~ 
0{\n.N/ ln(g — 1)), needs to be exponentially larger than any statistically relevant length 
scale {e.g. the coherence length). 

Another direction for further development would be to find a spin glass (or otherwise) 
model amenable to analytic treatment within the quantum cavity method. A "soft spin" 
Gaussian model would do, since the path integrals to be performed in an iteration could 
be computed exactly. However, we do not believe this model has a spin glass phase in the 
absence of higher order couplings. Whether one could treat such a coupling perturbatively 
is a question worthy of further exploration. 
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Note added: During the production of this manuscript, several other groups22i22 have 
independently developed quantum cavity techniques inspired by somewhat different ap- 



proaches to the quantum problem. Ref. 
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in several of these related techniques. In an important development, Ref. 



provides a quantitative comparison of the error 

describes a 
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continuous time version of the quantum cavity method described here which they apply to 



24 



a ferromagnetic model. 
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